function [U, S, V, timing] = blendenpik_sv(A, N, params)
% [x, timing] = blendenpik_over(A, n, params)
%
% Compute partial SVD decomposition: compute smallest n singular values
% of A. 
% 
% "params" - parameters governning the method.
%    params.type - type of mixing transform. Optional values: 'DCT', 'DHT', 'WHT'.
%                  Default is DHT.
%    params.gamma - gamma * n rows will be sampled (A is m-by-n). 
%                   Default is 4.
%    params.preprocess_steps - number of mixing steps to do in advance. 
%                               Default is 1.
%    params.maxcond - maximum condition number of the preconditioner.
%                     Default is 1 / (5 * epsilon_machine).
%    params.tol - convergence thershold for LOBPCG.
%                 Default is 1e-10.
%    params.maxit - maximum number of LOBPCG iterations.
%                   Default is 1000.
%
%
% Output:
%   U, S, V - as in SVD. 
%   timing - statistics on the time spent on various phases.
%          
% October 03, 2011


if (nargin < 3)
    params = struct;
end

if (~isfield(params, 'type'));
    params.type = 'DHT';
end

if (~isfield(params, 'gamma'))
    params.gamma = 4;
end

if (~isfield(params, 'preprocess_steps'))
    params.preprocess_steps = 1;
end

if (~isfield(params, 'tol'))
    params.tol = 1e-10;
end

if (~isfield(params, 'maxit'))
    params.maxit = 1000;
end

if (~isfield(params, 'maxcond'))
    params.maxcond = 1 / (5 * eps);
end


if (~isfield(params, 'slight_coherence'))
    params.slight_coherence = 0;
end

% TODO: improve start position using SVD of preconditioner. 

global U;
global F;

tstart = wtime;
[m, n] = size(A);

%% Build preconditioner
t1 = wtime;
[R, flag, timing.precond_timing, x0] = random_sample_precond(A, params);
timing.precond_total_time = wtime - t1;
disp(sprintf('\tBuilding preconditioner time: %.2f sec', timing.precond_total_time));

% TODO better use SVD instead of QR for rank-deficient and for more
% efficient gemm in iterations.

%% Solve
if (flag)  
    t1 = wtime;
    U = R;
    %U = eye(n);
    F = A;
    [V, S, flag, hs, hres] = lobpcg(rand(n, N), 'multF', [], 'precondU', [], ...
                                params.tol, params.maxit);
    if (flag ~= 0)
        % TODO
        disp(sprintf('\tLOBPCG failed!', timing.lsqr_time));
    else 
        S = sqrt(S);
        timing.lobpcg_time = wtime - t1;
        disp(sprintf('\tLOBPCG time: %.2f sec', timing.lobpcg_time));
        disp(sprintf('\tLOBPCG iterations: %d', size(hs, 2)));
    end
else
    disp(sprintf('Failed to get a full rank preconditioner. Using LAPACK.'));
    % TODO!
end

% TODO: get U
U = 0;

timing.total_time = wtime - tstart;
disp(sprintf('Total time: %.2f sec', timing.total_time));
